Behavioral Data

Demographics

First, let’s load the data and transform it into a Wide format table

behav <- read_tsv("data/behav/behav_data.txt",
                  col_types=cols())

behav$material[behav$material=="Swahili"] <- "Vocabulary"
behav$material[behav$material=="US Maps"] <- "Maps"

behav <- behav %>% rename(Gender = gender, Age = age)

wbehav <- behav %>% pivot_wider(values_from = alpha,
                                id_cols = c(Subject, Gender, Age),
                                names_from = material)

At this point, we can look at the participant demographocs

ggplot(data=wbehav, aes(x=Age, col=Gender)) +
  geom_histogram(aes(col="white", fill=Gender),
                 colour="white", alpha=0.5, 
                 position="identity", binwidth = 1) +
  scale_color_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  ggtitle("Age Distribution by Gender") + 
  theme_pander()
## Warning: Removed 2 rows containing non-finite values (stat_bin).

Rate of Forgetting (Alpha parameter) for the two tasks

We can now examine the distribution and correlation of the rate of forgetting for the verbal (Vocabulary) and Maps:

mu <- behav %>% 
  group_by(material) %>% 
  summarize(alpha=mean(alpha))
## `summarise()` ungrouping output (override with `.groups` argument)
p1 <- ggplot(behav, aes(x=alpha, fill=material)) +
  geom_histogram(col="white", binwidth = 0.025, 
                 alpha=0.5, position="identity") +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  geom_vline(data=mu, aes(xintercept=alpha, color=material),
             linetype="dashed") +
  xlab(expression(paste("Estimated value of ", alpha))) +
  theme_pander() +
  ylab("Count") +
  ggtitle("(A) Distribution By Materials") +
  theme(legend.position = c(0.2, 0.8)) +
  theme(plot.title = element_text(hjust = 0.5))

p2 <- ggplot(wbehav, aes(x=Vocabulary, y=Maps)) +
  geom_point(size=4, alpha=0.5, col="black") + #col=K[7]) +
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() + 
  scale_y_continuous() + 
  ggtitle("(B) Correlation Across Materials") +
  xlab(expression(paste(alpha, " Vocabulary"))) +
  ylab(expression(paste(alpha, " Maps"))) +
  geom_rug(col="black", lwd=1, alpha=.5) +
  annotate("segment", x=0.1, y=0.1, xend=0.5, 
           yend=0.5, col="grey", lwd=1, lty=2) +
  theme(plot.title = element_text(hjust = 0.5))
  
grid.arrange(p1, p2, ncol=2)

p3 <- ggplot(wbehav, aes(x=Vocabulary, y=Maps)) +
  geom_point(size=4, alpha=0.5, col="black") + #col=K[7]) +
  geom_smooth(method = "lm", formula = y ~ x, col="red", 
              fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() + 
  scale_y_continuous() + 
  coord_fixed() +
  ggtitle("Correlation of Alpha Parameters\nAcross Materials") +
  xlab(expression(paste(alpha, " Vocabulary"))) +
  ylab(expression(paste(alpha, " Maps"))) +
  annotate("segment", x=0.1, y=0.1, xend=0.5, 
           yend=0.5, col="grey", lwd=1, lty=2) +
  annotate("text", label = paste("R =", round(cor(wbehav$Vocabulary,
                                                  wbehav$Maps),
                                              2)), 
                                 x=0.2, y=0.45, col="red") +
  
  theme(plot.title = element_text(hjust = 0.5))

ggMarginal(p3, col="white", fill="black", alpha=0.5,
           type="hist", bins=13)

And here is a summary of the data:

summary(wbehav) %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Subject Gender Age Vocabulary Maps
X Length:50 Length:50 Min. :18.00 Min. :0.1219 Min. :0.2120
X.1 Class :character Class :character 1st Qu.:18.00 1st Qu.:0.2647 1st Qu.:0.3039
X.2 Mode :character Mode :character Median :19.00 Median :0.2917 Median :0.3447
X.3 NA NA Mean :18.92 Mean :0.2825 Mean :0.3355
X.4 NA NA 3rd Qu.:19.00 3rd Qu.:0.3203 3rd Qu.:0.3680
X.5 NA NA Max. :22.00 Max. :0.3852 Max. :0.4701
X.6 NA NA NA’s :2 NA NA

Eyes-Closed EEG Data

Now we look at the neural correlates of the rates of forgetting in EEG data. We begin with eyes closed EEG data, which is the most common type of resting-state EEG recording.

(Optional) Preprocess the Data

Set the process_raw variable to TRUE to reprocess the raw EEG data. That might take a very long time.

process_raw = F
if (process_raw) {
  setwd("./data/eyes_closed/")
  for (sub in dir()[grep("A[1-9]", dir())]) {
    setwd(sub)
    analyze.logfile(sub, "closed")
    setwd("..")
  }
  setwd("../..")
}

Load the Eyes-Closed Spectra and Summary Files

The preprocessing step produces a number of text files. Here, we are interested in the summary file, which contains many useful statistics, including the Individual Alpha Frequency (i.e., the Alpha Peak) for each channel, the spectra file, which contains the estimated log power for each channel in increments of 0.5 Hz.

ec_spectra <- NULL
for (sub in dir("data/eyes_closed/")[grep("A[1-9]",
                                        dir("data/eyes_closed/"))]) {
  table <- read_tsv(paste("data/eyes_closed/", sub, "/", 
                          sub, "_closed_spectra.txt", sep=""),
                    col_types = cols())
  if (is.null(ec_spectra)) {
    ec_spectra <- table
  } else {
    ec_spectra <- ec_spectra %>% bind_rows(table)
  }
}

ec_summary <- NULL
for (sub in dir("data/eyes_closed/")[grep("A[1-9]", 
                                          dir("data/eyes_closed/"))]) {
  table <- read_tsv(paste("data/eyes_closed/", sub, "/", 
                          sub, "_closed_summary.txt", sep=""),
                    col_types = cols())
  if (is.null(ec_summary)) {
    ec_summary <- table
  } else {
    ec_summary <- ec_summary %>% bind_rows(table)
  }
}

Individual Alpha frquency

The IAF is perhaps the most important subject level characteristic of the resting state spectrum. For this reason, we want to get some basic statistics. First, we calculate the subject-level IAF by computing the mode IAF across all channels for a given participant.

iaf_cols <- paste(c("AF3", "AF4", "F3", "F4", 
                    "F7", "F8", "FC5", "FC6", 
                    "T7", "T8", "P7", "P8", 
                    "O1", "O2"), 
                  "IAF", 
                  sep="_")

ec_iafs <- ec_summary %>%
  select(c("Subject", iaf_cols)) %>%
  filter(Subject %in% wbehav$Subject) %>%
  pivot_longer(cols = iaf_cols, names_sep="_", 
               names_to = c("Channel", "Discard"), 
               values_to = "IAF_EC") %>%
  select(Subject, Channel, IAF_EC) %>%
  group_by(Subject) %>%
  summarize(IAF_EC = Mode(IAF_EC))
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(iaf_cols)` instead of `iaf_cols` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## `summarise()` ungrouping output (override with `.groups` argument)

The mean IAF thus calculated is 10.3, which is exactly in the middle of our alpha frequency band.

ggplot(data=ec_iafs, aes(x = IAF_EC)) +
  geom_histogram(aes(col="white"), fill="black", 
                 colour="white", alpha=0.5, 
                 position="identity", binwidth = 0.5) +
  ggtitle("IAF Distribution, Eyes Closed") + 
  xlab("IAF (Eyes Closed)") +
  theme_pander()

Average Spectrum

First, let’s examine the spectrograms for each channel to make sure it looks normal and the IAF (individual alpha frequency) is reasonable, i.e., the same for each channel and roughly in the middle of the band definition

l_ec_spectra <- pivot_longer(ec_spectra, cols=names(ec_spectra)[3:130], 
                         names_to="Freq")

l_ec_spectra <- l_ec_spectra %>% 
  rename(Power = value) %>%
  add_column(Recording = "Eyes Closed")

l_ec_spectra$Frequency <- as.numeric(substr(l_ec_spectra$Freq, 
                                            0, str_length(l_ec_spectra$Freq) -2))
l_ec_spectra <- l_ec_spectra %>% 
  add_column(Band="Delta", BandMin=0, BandMax=4) 
l_ec_spectra$Band[l_ec_spectra$Frequency <= 40] <- "Gamma"
l_ec_spectra$Band[l_ec_spectra$Frequency < 30] <- "High Beta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 18] <- "Upper Beta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 15] <- "Low Beta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 13] <- "Alpha"
l_ec_spectra$Band[l_ec_spectra$Frequency < 8] <- "Theta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 4] <- "Delta"

l_ec_spectra$BandMin[l_ec_spectra$Frequency <= 40] <- 30
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 30] <- 18
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 18] <- 15
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 15] <- 13
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 13] <- 8
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 8] <- 4
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 4] <- 0

l_ec_spectra$BandMax[l_ec_spectra$Frequency <= 40] <- 40
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 30] <- 30
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 18] <- 18
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 15] <- 15
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 13] <- 13
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 8] <- 8
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 4] <- 4

l_ec_spectra$Band <- factor(l_ec_spectra$Band, 
                        levels = c("Delta", "Theta", "Alpha", 
                                   "Low Beta", "Upper Beta",
                                   "High Beta", "Gamma"))

Now, we remove the three participants for which we have poor quality data

l_ec_spectra <- l_ec_spectra %>%
  filter(Subject %in% behav$Subject)

We can visualize the power spectra to ensure that our data looks normal

gd <- l_ec_spectra %>% 
        group_by(Band) %>% 
        summarise(
          Min = mean(BandMin),
          Max = mean(BandMax),
          Power =mean(Power),
          Frequency = mean(BandMin)
          #Channel=
        )
## `summarise()` ungrouping output (override with `.groups` argument)
ggplot(data=l_ec_spectra, aes(x=Frequency, y=Power, Channel)) +
  geom_rect(data = gd, aes(xmin = Min, xmax = Max, fill = Band), 
            ymin=0, ymax=Inf, colour=NA, alpha=0.5) +
  stat_summary(fun.data=mean_sdl, 
               geom = "ribbon", colour = "white", alpha = 0.4) +
  stat_summary(fun = mean, geom = "line", lwd = 1) +
  facet_wrap(~ Channel, ncol=4) +
  scale_alpha_manual(values = seq(0.1, 0.9, 0.1)) +
  ggtitle("Eyes-Closed Power Spectrum Across Channels") +
  ylab("Log Power") + 
  xlab("Frequency (Hz)") + 
  theme_pander() +
  coord_cartesian(xlim=c(1,40), ylim=c(5, 17)) + #ylim=c(5,17)) +
  scale_fill_brewer(palette = "Set3") +
  theme(plot.title = element_text(hjust = 0.5))
## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

Now, let’s create a mean spectrogram with annotations that mark the different frequency bands. It serves as a visual reference that our data is not grossly out of whack.

K = brewer.pal(7, "Set3")

al_ec_spectra <- aggregate(l_ec_spectra[c("Power")],
                       list(Subject = l_ec_spectra$Subject, 
                            Frequency = l_ec_spectra$Frequency, 
                            Band = l_ec_spectra$Band), 
                       mean)

ggplot(data=al_ec_spectra, aes(x=Frequency, y=Power, col=Power)) +
  annotate("rect", xmin = 0, xmax = 4, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[1]) +
  annotate("rect", xmin = 4, xmax = 8, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[2]) +
  annotate("rect", xmin = 8, xmax = 13, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[3]) +      
  annotate("rect", xmin = 13, xmax = 15, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[4]) +
  annotate("rect", xmin = 15, xmax = 18, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill = K[5]) +
  annotate("rect", xmin = 18, xmax = 30, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[6]) +
  annotate("rect", xmin = 30, xmax = 40, ymin = 0, ymax = Inf,
           alpha = 0.5, fill=K[7]) +
  stat_summary(fun.data = mean_sdl, geom = "ribbon", 
               alpha = 0.5, fill = "grey50", col = "white") +
  stat_summary(fun=mean, geom="line", lwd=2) +
  coord_cartesian(xlim=c(0,40), ylim=c(5,18)) +
  xlab("Frequency (Hz)") +
  theme_pander() +
  annotate("text", x=2, y=5, label="Delta", angle=90, hjust=0) +
  annotate("text", x=6, y=5, label="Theta", angle=90, hjust=0) +
  annotate("text", x=10.5, y=5, label="Alpha", angle=90, hjust=0) +
  annotate("text", x=14, y=5, label="Low Beta", angle=90, hjust=0) +
  annotate("text", x=16.5, y=5, label="Upper Beta", angle=90, hjust=0) +
  annotate("text", x=24, y=5, label="High Beta", angle=90, hjust=0) +
  annotate("text", x=35, y=5, label="Gamma", angle=90, hjust=0) 
## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

Correlations with Rate of Forgetting

Having ensured that our EEG recordings look normal, reflect known neurophysiological characteristics, and that our frequency bands boundaries are reasonable, we can proceed with correlating the rates of forgetting with resting state QEEG characteristics.

First, we calculate the mean power for every band and channel, remove participants who were not included in the behavioral tests, and correct for multiple comparisons by frequency band:

wbehav <- wbehav %>%
  mutate(GlobalAlpha=(Vocabulary + Maps)/2)

Adata_ec <- l_ec_spectra %>% 
  group_by(Subject, Channel, Band, Recording) %>% 
  summarize(Power=mean(Power))
## `summarise()` regrouping output by 'Subject', 'Channel', 'Band' (override with `.groups` argument)
Fdata_ec <- inner_join(Adata_ec, behav, by="Subject")

Rdata_ec <- Fdata_ec %>%
  group_by(material, Channel, Band, Recording) %>%
  summarise(r = cor(Power, alpha),
            p = cor.test(Power, alpha)$p.value)
## `summarise()` regrouping output by 'material', 'Channel', 'Band' (override with `.groups` argument)
Rdata_ec <- Rdata_ec %>%
  group_by(material, Band, Recording) %>%
  mutate(q = p.adjust(p, method="BH")) 

Here is the full table of statistics:

write_csv(Rdata_ec, "correlations_channel_band_eyes_closed.csv", col_names = T)

Rdata_ec %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
material Channel Band Recording r p q
Maps AF3 Delta Eyes Closed 0.010 0.947 0.971
Maps AF3 Theta Eyes Closed 0.072 0.621 0.790
Maps AF3 Alpha Eyes Closed 0.113 0.436 0.592
Maps AF3 Low Beta Eyes Closed 0.189 0.188 0.316
Maps AF3 Upper Beta Eyes Closed 0.314 0.026 0.092
Maps AF3 High Beta Eyes Closed 0.219 0.127 0.403
Maps AF3 Gamma Eyes Closed 0.101 0.484 0.924
Maps AF4 Delta Eyes Closed 0.096 0.507 0.971
Maps AF4 Theta Eyes Closed 0.079 0.585 0.790
Maps AF4 Alpha Eyes Closed 0.147 0.308 0.582
Maps AF4 Low Beta Eyes Closed 0.216 0.132 0.316
Maps AF4 Upper Beta Eyes Closed 0.339 0.016 0.092
Maps AF4 High Beta Eyes Closed 0.280 0.049 0.342
Maps AF4 Gamma Eyes Closed 0.145 0.315 0.924
Maps F3 Delta Eyes Closed -0.005 0.971 0.971
Maps F3 Theta Eyes Closed 0.072 0.619 0.790
Maps F3 Alpha Eyes Closed 0.140 0.333 0.582
Maps F3 Low Beta Eyes Closed 0.183 0.203 0.316
Maps F3 Upper Beta Eyes Closed 0.256 0.073 0.120
Maps F3 High Beta Eyes Closed 0.184 0.200 0.403
Maps F3 Gamma Eyes Closed 0.022 0.881 0.924
Maps F4 Delta Eyes Closed 0.239 0.095 0.343
Maps F4 Theta Eyes Closed 0.122 0.400 0.790
Maps F4 Alpha Eyes Closed 0.156 0.278 0.582
Maps F4 Low Beta Eyes Closed 0.200 0.164 0.316
Maps F4 Upper Beta Eyes Closed 0.290 0.041 0.116
Maps F4 High Beta Eyes Closed 0.215 0.133 0.403
Maps F4 Gamma Eyes Closed 0.135 0.349 0.924
Maps F7 Delta Eyes Closed -0.058 0.689 0.971
Maps F7 Theta Eyes Closed 0.087 0.550 0.790
Maps F7 Alpha Eyes Closed 0.115 0.425 0.592
Maps F7 Low Beta Eyes Closed 0.203 0.158 0.316
Maps F7 Upper Beta Eyes Closed 0.254 0.075 0.120
Maps F7 High Beta Eyes Closed 0.121 0.403 0.564
Maps F7 Gamma Eyes Closed -0.021 0.885 0.924
Maps F8 Delta Eyes Closed 0.011 0.939 0.971
Maps F8 Theta Eyes Closed 0.147 0.309 0.790
Maps F8 Alpha Eyes Closed 0.173 0.230 0.582
Maps F8 Low Beta Eyes Closed 0.206 0.152 0.316
Maps F8 Upper Beta Eyes Closed 0.319 0.024 0.092
Maps F8 High Beta Eyes Closed 0.184 0.201 0.403
Maps F8 Gamma Eyes Closed 0.065 0.655 0.924
Maps FC5 Delta Eyes Closed -0.030 0.838 0.971
Maps FC5 Theta Eyes Closed 0.025 0.865 0.932
Maps FC5 Alpha Eyes Closed 0.078 0.591 0.637
Maps FC5 Low Beta Eyes Closed 0.128 0.376 0.526
Maps FC5 Upper Beta Eyes Closed 0.153 0.288 0.367
Maps FC5 High Beta Eyes Closed 0.071 0.622 0.725
Maps FC5 Gamma Eyes Closed 0.014 0.924 0.924
Maps FC6 Delta Eyes Closed 0.069 0.633 0.971
Maps FC6 Theta Eyes Closed 0.106 0.462 0.790
Maps FC6 Alpha Eyes Closed 0.167 0.247 0.582
Maps FC6 Low Beta Eyes Closed 0.183 0.203 0.316
Maps FC6 Upper Beta Eyes Closed 0.272 0.056 0.120
Maps FC6 High Beta Eyes Closed 0.153 0.290 0.507
Maps FC6 Gamma Eyes Closed 0.088 0.545 0.924
Maps O1 Delta Eyes Closed 0.221 0.123 0.343
Maps O1 Theta Eyes Closed 0.009 0.952 0.952
Maps O1 Alpha Eyes Closed -0.096 0.508 0.592
Maps O1 Low Beta Eyes Closed -0.039 0.786 0.847
Maps O1 Upper Beta Eyes Closed 0.100 0.491 0.573
Maps O1 High Beta Eyes Closed 0.093 0.519 0.660
Maps O1 Gamma Eyes Closed 0.139 0.336 0.924
Maps O2 Delta Eyes Closed 0.229 0.109 0.343
Maps O2 Theta Eyes Closed 0.106 0.463 0.790
Maps O2 Alpha Eyes Closed 0.060 0.680 0.680
Maps O2 Low Beta Eyes Closed 0.098 0.497 0.632
Maps O2 Upper Beta Eyes Closed 0.193 0.180 0.252
Maps O2 High Beta Eyes Closed 0.205 0.153 0.403
Maps O2 Gamma Eyes Closed 0.181 0.209 0.924
Maps P7 Delta Eyes Closed 0.201 0.161 0.375
Maps P7 Theta Eyes Closed -0.037 0.798 0.931
Maps P7 Alpha Eyes Closed -0.148 0.306 0.582
Maps P7 Low Beta Eyes Closed -0.028 0.847 0.847
Maps P7 Upper Beta Eyes Closed 0.049 0.735 0.791
Maps P7 High Beta Eyes Closed 0.027 0.852 0.852
Maps P7 Gamma Eyes Closed 0.028 0.846 0.924
Maps P8 Delta Eyes Closed 0.326 0.021 0.290
Maps P8 Theta Eyes Closed 0.343 0.015 0.208
Maps P8 Alpha Eyes Closed 0.222 0.122 0.582
Maps P8 Low Beta Eyes Closed 0.298 0.036 0.316
Maps P8 Upper Beta Eyes Closed 0.365 0.009 0.092
Maps P8 High Beta Eyes Closed 0.332 0.018 0.258
Maps P8 Gamma Eyes Closed 0.288 0.043 0.599
Maps T7 Delta Eyes Closed -0.042 0.770 0.971
Maps T7 Theta Eyes Closed -0.148 0.307 0.790
Maps T7 Alpha Eyes Closed -0.099 0.494 0.592
Maps T7 Low Beta Eyes Closed -0.060 0.681 0.795
Maps T7 Upper Beta Eyes Closed -0.023 0.876 0.876
Maps T7 High Beta Eyes Closed -0.037 0.797 0.852
Maps T7 Gamma Eyes Closed -0.017 0.905 0.924
Maps T8 Delta Eyes Closed 0.247 0.084 0.343
Maps T8 Theta Eyes Closed 0.251 0.078 0.547
Maps T8 Alpha Eyes Closed 0.239 0.095 0.582
Maps T8 Low Beta Eyes Closed 0.214 0.135 0.316
Maps T8 Upper Beta Eyes Closed 0.252 0.077 0.120
Maps T8 High Beta Eyes Closed 0.140 0.332 0.517
Maps T8 Gamma Eyes Closed 0.036 0.804 0.924
Vocabulary AF3 Delta Eyes Closed -0.008 0.958 0.958
Vocabulary AF3 Theta Eyes Closed 0.023 0.873 0.962
Vocabulary AF3 Alpha Eyes Closed 0.086 0.553 0.779
Vocabulary AF3 Low Beta Eyes Closed 0.298 0.035 0.145
Vocabulary AF3 Upper Beta Eyes Closed 0.225 0.116 0.334
Vocabulary AF3 High Beta Eyes Closed 0.020 0.890 0.958
Vocabulary AF3 Gamma Eyes Closed -0.091 0.529 0.617
Vocabulary AF4 Delta Eyes Closed -0.059 0.683 0.939
Vocabulary AF4 Theta Eyes Closed 0.036 0.802 0.962
Vocabulary AF4 Alpha Eyes Closed 0.105 0.466 0.779
Vocabulary AF4 Low Beta Eyes Closed 0.257 0.071 0.157
Vocabulary AF4 Upper Beta Eyes Closed 0.223 0.119 0.334
Vocabulary AF4 High Beta Eyes Closed 0.026 0.859 0.958
Vocabulary AF4 Gamma Eyes Closed -0.114 0.432 0.605
Vocabulary F3 Delta Eyes Closed -0.201 0.162 0.939
Vocabulary F3 Theta Eyes Closed -0.007 0.962 0.962
Vocabulary F3 Alpha Eyes Closed 0.064 0.660 0.840
Vocabulary F3 Low Beta Eyes Closed 0.235 0.101 0.157
Vocabulary F3 Upper Beta Eyes Closed 0.175 0.223 0.428
Vocabulary F3 High Beta Eyes Closed -0.035 0.808 0.958
Vocabulary F3 Gamma Eyes Closed -0.191 0.185 0.581
Vocabulary F4 Delta Eyes Closed -0.041 0.778 0.939
Vocabulary F4 Theta Eyes Closed 0.036 0.802 0.962
Vocabulary F4 Alpha Eyes Closed 0.105 0.469 0.779
Vocabulary F4 Low Beta Eyes Closed 0.273 0.055 0.154
Vocabulary F4 Upper Beta Eyes Closed 0.249 0.082 0.334
Vocabulary F4 High Beta Eyes Closed -0.001 0.997 0.997
Vocabulary F4 Gamma Eyes Closed -0.140 0.332 0.581
Vocabulary F7 Delta Eyes Closed 0.015 0.917 0.958
Vocabulary F7 Theta Eyes Closed 0.029 0.842 0.962
Vocabulary F7 Alpha Eyes Closed 0.099 0.494 0.779
Vocabulary F7 Low Beta Eyes Closed 0.290 0.041 0.145
Vocabulary F7 Upper Beta Eyes Closed 0.168 0.245 0.428
Vocabulary F7 High Beta Eyes Closed -0.043 0.765 0.958
Vocabulary F7 Gamma Eyes Closed -0.126 0.383 0.595
Vocabulary F8 Delta Eyes Closed -0.126 0.383 0.939
Vocabulary F8 Theta Eyes Closed -0.019 0.897 0.962
Vocabulary F8 Alpha Eyes Closed 0.137 0.343 0.779
Vocabulary F8 Low Beta Eyes Closed 0.241 0.091 0.157
Vocabulary F8 Upper Beta Eyes Closed 0.196 0.173 0.403
Vocabulary F8 High Beta Eyes Closed -0.082 0.571 0.958
Vocabulary F8 Gamma Eyes Closed -0.206 0.151 0.581
Vocabulary FC5 Delta Eyes Closed -0.087 0.546 0.939
Vocabulary FC5 Theta Eyes Closed -0.090 0.532 0.962
Vocabulary FC5 Alpha Eyes Closed -0.007 0.962 0.962
Vocabulary FC5 Low Beta Eyes Closed 0.148 0.306 0.357
Vocabulary FC5 Upper Beta Eyes Closed 0.046 0.751 0.809
Vocabulary FC5 High Beta Eyes Closed -0.115 0.428 0.958
Vocabulary FC5 Gamma Eyes Closed -0.150 0.297 0.581
Vocabulary FC6 Delta Eyes Closed -0.092 0.524 0.939
Vocabulary FC6 Theta Eyes Closed -0.040 0.781 0.962
Vocabulary FC6 Alpha Eyes Closed 0.097 0.501 0.779
Vocabulary FC6 Low Beta Eyes Closed 0.203 0.157 0.214
Vocabulary FC6 Upper Beta Eyes Closed 0.154 0.286 0.444
Vocabulary FC6 High Beta Eyes Closed -0.084 0.562 0.958
Vocabulary FC6 Gamma Eyes Closed -0.172 0.232 0.581
Vocabulary O1 Delta Eyes Closed -0.044 0.759 0.939
Vocabulary O1 Theta Eyes Closed 0.057 0.697 0.962
Vocabulary O1 Alpha Eyes Closed 0.041 0.776 0.866
Vocabulary O1 Low Beta Eyes Closed 0.198 0.168 0.214
Vocabulary O1 Upper Beta Eyes Closed 0.138 0.340 0.476
Vocabulary O1 High Beta Eyes Closed -0.052 0.722 0.958
Vocabulary O1 Gamma Eyes Closed -0.145 0.315 0.581
Vocabulary O2 Delta Eyes Closed -0.036 0.804 0.939
Vocabulary O2 Theta Eyes Closed 0.077 0.597 0.962
Vocabulary O2 Alpha Eyes Closed 0.122 0.399 0.779
Vocabulary O2 Low Beta Eyes Closed 0.340 0.016 0.111
Vocabulary O2 Upper Beta Eyes Closed 0.271 0.057 0.334
Vocabulary O2 High Beta Eyes Closed 0.071 0.625 0.958
Vocabulary O2 Gamma Eyes Closed -0.092 0.524 0.617
Vocabulary P7 Delta Eyes Closed -0.065 0.656 0.939
Vocabulary P7 Theta Eyes Closed -0.121 0.404 0.962
Vocabulary P7 Alpha Eyes Closed -0.085 0.556 0.779
Vocabulary P7 Low Beta Eyes Closed 0.040 0.783 0.843
Vocabulary P7 Upper Beta Eyes Closed -0.020 0.892 0.892
Vocabulary P7 High Beta Eyes Closed -0.119 0.410 0.958
Vocabulary P7 Gamma Eyes Closed -0.184 0.201 0.581
Vocabulary P8 Delta Eyes Closed 0.072 0.620 0.939
Vocabulary P8 Theta Eyes Closed 0.226 0.115 0.962
Vocabulary P8 Alpha Eyes Closed 0.210 0.142 0.779
Vocabulary P8 Low Beta Eyes Closed 0.424 0.002 0.030
Vocabulary P8 Upper Beta Eyes Closed 0.335 0.017 0.242
Vocabulary P8 High Beta Eyes Closed 0.158 0.273 0.958
Vocabulary P8 Gamma Eyes Closed 0.013 0.926 0.926
Vocabulary T7 Delta Eyes Closed -0.051 0.726 0.939
Vocabulary T7 Theta Eyes Closed -0.104 0.471 0.962
Vocabulary T7 Alpha Eyes Closed 0.036 0.804 0.866
Vocabulary T7 Low Beta Eyes Closed -0.027 0.851 0.851
Vocabulary T7 Upper Beta Eyes Closed -0.059 0.684 0.798
Vocabulary T7 High Beta Eyes Closed -0.058 0.690 0.958
Vocabulary T7 Gamma Eyes Closed -0.064 0.660 0.710
Vocabulary T8 Delta Eyes Closed -0.037 0.798 0.939
Vocabulary T8 Theta Eyes Closed 0.041 0.779 0.962
Vocabulary T8 Alpha Eyes Closed 0.173 0.230 0.779
Vocabulary T8 Low Beta Eyes Closed 0.238 0.096 0.157
Vocabulary T8 Upper Beta Eyes Closed 0.097 0.504 0.641
Vocabulary T8 High Beta Eyes Closed -0.096 0.506 0.958
Vocabulary T8 Gamma Eyes Closed -0.225 0.116 0.581

And here is the corresponding distributions of r and p values. The dashed lines correspond to a significant threshold of p < .05 on either a r or a p-value scale.

ggplot(Rdata_ec, aes(x = Band, y = r, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", 
               alpha=0.5, geom = "errorbar") +
  facet_wrap(~ material) +
  ggtitle("Correlation by Frequency Band") +
  ylab("r value") +
  annotate("segment", x=-Inf, xend=Inf, y=0.28, yend = 0.28, lty=2) +
  annotate("segment", x=-Inf, xend=Inf, y=-0.28, yend = -0.28, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  theme_pander()

ggplot(Rdata_ec, aes(x = Band, y = p, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", alpha=0.5, geom = "errorbar") +
  ggtitle("p-value, by Band") +
  ylab("p value") +
  facet_wrap(~ material) +
  annotate("segment", x=-Inf, xend=Inf, y=0.05, yend = 0.05, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  scale_y_log10() +
  theme_pander()

From the original list, we can extract only those channels that survive the FDR correction:

survivors <- Rdata_ec %>%
  filter(q < 0.05)
  
survivors %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
material Channel Band Recording r p q
Vocabulary P8 Low Beta Eyes Closed 0.424 0.002 0.03

Only one channel P8, in one frequency band (low beta) survives FDR correction. We can plot the correlation between power in P8 and rate of forgetting to see the relationship:

focus <- Fdata_ec %>%
  filter(material == "Vocabulary", 
         Channel %in% survivors$Channel,
         Band == "Low Beta") %>%
  rename(Alpha = alpha)

p <- ggplot(focus, aes(x = Alpha, y = Power)) +
  geom_point(size = 4, alpha = 0.5, col = "black") + 
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() +   
  scale_y_continuous() + 
  ggtitle("Rate of Forgetting\nand Eyes-Closed Beta Power") +
  xlab(expression(paste(alpha, " Vocabulary"))) +
  ylab("Low Beta (13-15 Hz) Power Over P8") +
  geom_text(data=survivors, col="red",
            mapping=aes(x=0.25, y=10.5, 
                        label= paste("r =", round(r, 2)))) +
  theme(plot.title = element_text(hjust = 0.5))

ggMarginal(p, col = "white", fill = "black", alpha = 0.5,
           type = "hist", bins = 13)

Eyes-Open EEG Data

(Optional) Preprocess the Data

Set the process_raw variable to TRUE to reprocess the raw EEG data. That might take a very long time.

process_raw = F
if (process_raw) {
  setwd("./data/eyes_open/")
  for (sub in dir()[grep("A[1-9]", dir())]) {
    setwd(sub)
    analyze.logfile(sub, "open")
    setwd("..")
  }
  setwd("../..")
}

Load the Eyes-Open Spectra and Summary Files

As in the case of eyes-closed data, we are going to examine the spectral and IAF characteristics of eyes-open recordings though the summary files produced by preprocessing.

eo_spectra <- NULL
for (sub in dir("data/eyes_open/")[grep("A[1-9]", 
                                        dir("data/eyes_open/"))]) {
  table <- read_tsv(paste("data/eyes_open/", sub, "/", 
                          sub, "_open_spectra.txt", sep=""),
                    col_types = cols())
  if (is.null(eo_spectra)) {
    eo_spectra <- table
  } else {
    eo_spectra <- eo_spectra %>% bind_rows(table)
  }
}

eo_summary <- NULL
for (sub in dir("data/eyes_open/")[grep("A[1-9]", 
                                        dir("data/eyes_open/"))]) {
  table <- read_tsv(paste("data/eyes_open/", sub, "/", 
                          sub, "_open_summary.txt", sep=""),
                    col_types = cols())
  if (is.null(eo_summary)) {
    eo_summary <- table
  } else {
    eo_summary <- eo_summary %>% bind_rows(table)
  }
}

Individual Alpha frquency

Since eyes-open recordings are somewhat less common, we need to run a few sanity checks on them. The first and most obsvious concerns the IAF, and whether the IAFs during eyes-open recordings are similar and correlated to those observed during eyes-closed sessions.

eo_iafs <- eo_summary %>%
  select(c("Subject", iaf_cols)) %>%
  filter(Subject %in% wbehav$Subject) %>%
  pivot_longer(cols = iaf_cols, names_sep="_", 
               names_to = c("Channel", "Discard"), 
               values_to = "IAF_EO") %>%
  select(Subject, Channel, IAF_EO) %>%
  group_by(Subject) %>%
  summarize(IAF_EO = Mode(IAF_EO))
## `summarise()` ungrouping output (override with `.groups` argument)
iafs <- inner_join(eo_iafs, ec_iafs, by="Subject")

The mean IAF thus calculated is 10.33, which is, once more, exactly in the middle of our alpha frequency band.

ggplot(data=eo_iafs, aes(x = IAF_EO)) +
  geom_histogram(aes(col="white"), fill="black", 
                 colour="white", alpha=0.5, 
                 position="identity", binwidth = 0.5) +
  ggtitle("IAF Distribution, Eyes-Open") + 
  xlab("IAF (Eyes Open)") +
  theme_pander()

As it can be seen, there is a significant correlation in the IAFs between the two recordings (p = 10^{-6}). Between recordings, the individual IAF remains the same +/- 0.83.

ggplot(iafs, aes(x=IAF_EO, y=IAF_EC)) +
  geom_count(alpha=0.5, col="black") + 
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  xlab("IAF Eyes Open (Hz)") +
  ylab("IAF Eyes Closed (Hz") +
  geom_text(data = summarize(iafs, Mean=mean(IAF_EC)),  # Avoid label overlaps
            x= 9.5, y=11.5, 
            col="red", 
            label=paste("r =", 
                        round(cor(iafs$IAF_EC, iafs$IAF_EO), 2))) +
  ggtitle("Correlation Between Eyes-Closed and Eyes-Open IAF") +
  theme_pander()

Average Spectrum

We can also investigate whether the average spectrum is comparable to that of eyes-closed recordings. To do so, we will calculate the average spectrum across channels, much in the same ways as it was done for eyes-closed data:

l_eo_spectra <- pivot_longer(eo_spectra, cols=names(ec_spectra)[3:130], 
                         names_to="Freq")

l_eo_spectra <- l_eo_spectra %>% 
  rename(Power = value) %>%
  add_column(Recording = "Eyes Open")

l_eo_spectra$Frequency <- as.numeric(substr(l_eo_spectra$Freq, 
                                            0, str_length(l_eo_spectra$Freq) -2))
l_eo_spectra <- l_eo_spectra %>% 
  add_column(Band="Delta", BandMin=0, BandMax=4) 
l_eo_spectra$Band[l_eo_spectra$Frequency <= 40] <- "Gamma"
l_eo_spectra$Band[l_eo_spectra$Frequency < 30] <- "High Beta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 18] <- "Upper Beta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 15] <- "Low Beta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 13] <- "Alpha"
l_eo_spectra$Band[l_eo_spectra$Frequency < 8] <- "Theta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 4] <- "Delta"

l_eo_spectra$BandMin[l_eo_spectra$Frequency <= 40] <- 30
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 30] <- 18
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 18] <- 15
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 15] <- 13
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 13] <- 8
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 8] <- 4
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 4] <- 0

l_eo_spectra$BandMax[l_eo_spectra$Frequency <= 40] <- 40
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 30] <- 30
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 18] <- 18
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 15] <- 15
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 13] <- 13
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 8] <- 8
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 4] <- 4

l_eo_spectra$Band <- factor(l_eo_spectra$Band, 
                        levels = c("Delta", "Theta", "Alpha", 
                                   "Low Beta", "Upper Beta",
                                   "High Beta", "Gamma"))

Now, we remove the three participants for which we have poor quality data

l_eo_spectra <- l_eo_spectra %>%
  filter(Subject %in% behav$Subject)

And we can visualize the eyes-open power spectra by channel:

gd <- l_eo_spectra %>% 
        group_by(Band) %>% 
        summarise(
          Min = mean(BandMin),
          Max = mean(BandMax),
          Power =mean(Power),
          Frequency = mean(BandMin)
        )
## `summarise()` ungrouping output (override with `.groups` argument)
ggplot(data=l_eo_spectra, aes(x=Frequency, y=Power, Channel)) +
  geom_rect(data = gd, aes(xmin = Min, xmax = Max, fill = Band), 
            ymin=0, ymax=Inf, colour=NA, alpha=0.5) +
  stat_summary(fun.data=mean_sdl, 
               geom = "ribbon", colour = "white", alpha = 0.4) +
  stat_summary(fun = mean, geom = "line", lwd = 1) +
  facet_wrap(~ Channel, ncol=4) +
  scale_alpha_manual(values = seq(0.1, 0.9, 0.1)) +
  ggtitle("Eyes-Open Power Spectrum Across Channels") +
  ylab("Log Power") + 
  xlab("Frequency (Hz)") + 
  theme_pander() +
  coord_cartesian(xlim=c(1,40), ylim=c(5, 17)) + #ylim=c(5,17)) +
  scale_fill_brewer(palette = "Set3") +
  theme(plot.title = element_text(hjust = 0.5))
## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

At a first sight, the eyes-open power spectrum looks remarkably similar, with the obvious (and expected) difference of a less prominent pean in the alpha band. To better compare the two types of recordings, we can visualize the avearage spectra on top of ech other:

l_spectra <- rbind(l_ec_spectra, l_eo_spectra)

al_spectra <- l_spectra %>%
  group_by(Subject, Frequency, Band, Recording) %>%
  summarise(Power = mean(Power))
## `summarise()` regrouping output by 'Subject', 'Frequency', 'Band' (override with `.groups` argument)
al_spectra <- al_spectra %>% ungroup()

ggplot(data=al_spectra, aes(x = Frequency, y = Power, 
                            col = Recording, fill = Recording)) +
  annotate("rect", xmin = 0, xmax = 4, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[1]) +
  annotate("rect", xmin = 4, xmax = 8, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[2]) +
  annotate("rect", xmin = 8, xmax = 13, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[3]) +      
  annotate("rect", xmin = 13, xmax = 15, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[4]) +
  annotate("rect", xmin = 15, xmax = 18, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill = K[5]) +
  annotate("rect", xmin = 18, xmax = 30, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[6]) +
  annotate("rect", xmin = 30, xmax = 40, ymin = 0, ymax = Inf,
           alpha = 0.5, fill=K[7]) +
  stat_summary(fun.data = mean_sdl, geom = "ribbon", 
               alpha = 0.5, col = "white") +
  scale_color_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  stat_summary(fun = mean, geom = "line", lwd = 2) +
  coord_cartesian(xlim=c(0,40), ylim=c(5,18)) +
  xlab("Frequency (Hz)") +
  ylab("Log Power") +
  theme_pander() +
  annotate("text", x=2, y=18, label="Delta", angle=90, hjust=1) +
  annotate("text", x=6, y=18, label="Theta", angle=90, hjust=1) +
  annotate("text", x=10.5, y=18, label="Alpha", angle=90, hjust=1) +
  annotate("text", x=14, y=18, label="Low Beta", angle=90, hjust=1) +
  annotate("text", x=16.5, y=18, label="Upper Beta", angle=90, hjust=1) +
  annotate("text", x=24, y=18, label="High Beta", angle=90, hjust=1) +
  annotate("text", x=35, y=18, label="Gamma", angle=90, hjust=1) 
## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

Correlations with Rate of Forgetting

Like we did for eyes-closed data, we can now calculate the correlations between eyes-open EEG power spectra and the behavioral rate fo forgetting for verbal and visual materials.

wbehav <- wbehav %>%
  mutate(GlobalAlpha=(Vocabulary + Maps) / 2)

Adata_eo <- l_eo_spectra %>% 
  group_by(Subject, Channel, Band, Recording) %>% 
  summarize(Power=mean(Power))
## `summarise()` regrouping output by 'Subject', 'Channel', 'Band' (override with `.groups` argument)
Fdata_eo <- inner_join(Adata_eo, behav, by = "Subject")

Cdata_eo <- Fdata_eo %>%
  group_by(material, Channel, Band, Recording) %>%
  summarise(r = cor(Power, alpha),
            p = cor.test(Power, alpha)$p.value)
## `summarise()` regrouping output by 'material', 'Channel', 'Band' (override with `.groups` argument)
Rdata_eo <- Cdata_eo %>%
  group_by(material, Band, Recording) %>%
  mutate(q = p.adjust(p, method="BH")) 

We save the data so that it can be visualized as a topomap with Python MNE. Here is the full table of statistics:

write_csv(Rdata_eo, "correlations_channel_band_eyes_open.csv", col_names = T)

Rdata_eo %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
material Channel Band Recording r p q
Maps AF3 Delta Eyes Open 0.180 0.211 0.686
Maps AF3 Theta Eyes Open 0.122 0.397 0.556
Maps AF3 Alpha Eyes Open 0.218 0.129 0.200
Maps AF3 Low Beta Eyes Open 0.255 0.074 0.206
Maps AF3 Upper Beta Eyes Open 0.383 0.006 0.021
Maps AF3 High Beta Eyes Open 0.279 0.050 0.140
Maps AF3 Gamma Eyes Open 0.236 0.100 0.447
Maps AF4 Delta Eyes Open 0.141 0.330 0.686
Maps AF4 Theta Eyes Open 0.136 0.347 0.540
Maps AF4 Alpha Eyes Open 0.247 0.083 0.184
Maps AF4 Low Beta Eyes Open 0.301 0.034 0.206
Maps AF4 Upper Beta Eyes Open 0.383 0.006 0.021
Maps AF4 High Beta Eyes Open 0.315 0.026 0.140
Maps AF4 Gamma Eyes Open 0.256 0.073 0.447
Maps F3 Delta Eyes Open 0.065 0.656 0.707
Maps F3 Theta Eyes Open 0.246 0.085 0.441
Maps F3 Alpha Eyes Open 0.232 0.105 0.184
Maps F3 Low Beta Eyes Open 0.221 0.122 0.244
Maps F3 Upper Beta Eyes Open 0.326 0.021 0.037
Maps F3 High Beta Eyes Open 0.219 0.126 0.253
Maps F3 Gamma Eyes Open 0.126 0.384 0.662
Maps F4 Delta Eyes Open 0.103 0.478 0.686
Maps F4 Theta Eyes Open 0.236 0.099 0.441
Maps F4 Alpha Eyes Open 0.248 0.082 0.184
Maps F4 Low Beta Eyes Open 0.265 0.063 0.206
Maps F4 Upper Beta Eyes Open 0.338 0.016 0.034
Maps F4 High Beta Eyes Open 0.294 0.038 0.140
Maps F4 Gamma Eyes Open 0.204 0.154 0.447
Maps F7 Delta Eyes Open 0.233 0.104 0.686
Maps F7 Theta Eyes Open 0.148 0.304 0.540
Maps F7 Alpha Eyes Open 0.247 0.084 0.184
Maps F7 Low Beta Eyes Open 0.281 0.048 0.206
Maps F7 Upper Beta Eyes Open 0.428 0.002 0.013
Maps F7 High Beta Eyes Open 0.279 0.050 0.140
Maps F7 Gamma Eyes Open 0.202 0.160 0.447
Maps F8 Delta Eyes Open 0.071 0.624 0.707
Maps F8 Theta Eyes Open 0.175 0.225 0.526
Maps F8 Alpha Eyes Open 0.287 0.043 0.184
Maps F8 Low Beta Eyes Open 0.242 0.091 0.212
Maps F8 Upper Beta Eyes Open 0.336 0.017 0.034
Maps F8 High Beta Eyes Open 0.199 0.166 0.291
Maps F8 Gamma Eyes Open 0.088 0.544 0.662
Maps FC5 Delta Eyes Open 0.032 0.826 0.826
Maps FC5 Theta Eyes Open -0.028 0.845 0.985
Maps FC5 Alpha Eyes Open 0.167 0.247 0.315
Maps FC5 Low Beta Eyes Open 0.189 0.188 0.263
Maps FC5 Upper Beta Eyes Open 0.229 0.109 0.153
Maps FC5 High Beta Eyes Open 0.137 0.343 0.403
Maps FC5 Gamma Eyes Open 0.094 0.515 0.662
Maps FC6 Delta Eyes Open -0.115 0.425 0.686
Maps FC6 Theta Eyes Open 0.094 0.515 0.655
Maps FC6 Alpha Eyes Open 0.190 0.187 0.261
Maps FC6 Low Beta Eyes Open 0.132 0.362 0.461
Maps FC6 Upper Beta Eyes Open 0.179 0.214 0.249
Maps FC6 High Beta Eyes Open 0.136 0.346 0.403
Maps FC6 Gamma Eyes Open 0.037 0.798 0.798
Maps O1 Delta Eyes Open 0.074 0.608 0.707
Maps O1 Theta Eyes Open 0.011 0.938 0.988
Maps O1 Alpha Eyes Open 0.085 0.556 0.649
Maps O1 Low Beta Eyes Open 0.026 0.857 0.857
Maps O1 Upper Beta Eyes Open 0.191 0.183 0.233
Maps O1 High Beta Eyes Open 0.145 0.314 0.403
Maps O1 Gamma Eyes Open 0.126 0.383 0.662
Maps O2 Delta Eyes Open 0.123 0.395 0.686
Maps O2 Theta Eyes Open 0.137 0.342 0.540
Maps O2 Alpha Eyes Open 0.238 0.096 0.184
Maps O2 Low Beta Eyes Open 0.210 0.144 0.244
Maps O2 Upper Beta Eyes Open 0.344 0.015 0.034
Maps O2 High Beta Eyes Open 0.239 0.094 0.220
Maps O2 Gamma Eyes Open 0.176 0.222 0.518
Maps P7 Delta Eyes Open 0.128 0.374 0.686
Maps P7 Theta Eyes Open -0.002 0.988 0.988
Maps P7 Alpha Eyes Open 0.068 0.640 0.689
Maps P7 Low Beta Eyes Open 0.041 0.778 0.838
Maps P7 Upper Beta Eyes Open 0.112 0.440 0.474
Maps P7 High Beta Eyes Open 0.074 0.608 0.608
Maps P7 Gamma Eyes Open 0.062 0.669 0.720
Maps P8 Delta Eyes Open 0.117 0.419 0.686
Maps P8 Theta Eyes Open 0.275 0.053 0.441
Maps P8 Alpha Eyes Open 0.316 0.025 0.184
Maps P8 Low Beta Eyes Open 0.320 0.024 0.206
Maps P8 Upper Beta Eyes Open 0.436 0.002 0.013
Maps P8 High Beta Eyes Open 0.311 0.028 0.140
Maps P8 Gamma Eyes Open 0.237 0.098 0.447
Maps T7 Delta Eyes Open -0.100 0.490 0.686
Maps T7 Theta Eyes Open -0.191 0.184 0.515
Maps T7 Alpha Eyes Open -0.033 0.818 0.818
Maps T7 Low Beta Eyes Open -0.092 0.524 0.611
Maps T7 Upper Beta Eyes Open -0.072 0.621 0.621
Maps T7 High Beta Eyes Open -0.137 0.344 0.403
Maps T7 Gamma Eyes Open -0.105 0.468 0.662
Maps T8 Delta Eyes Open 0.198 0.168 0.686
Maps T8 Theta Eyes Open 0.219 0.126 0.441
Maps T8 Alpha Eyes Open 0.270 0.058 0.184
Maps T8 Low Beta Eyes Open 0.203 0.157 0.244
Maps T8 Upper Beta Eyes Open 0.234 0.102 0.153
Maps T8 High Beta Eyes Open 0.081 0.576 0.608
Maps T8 Gamma Eyes Open 0.083 0.567 0.662
Vocabulary AF3 Delta Eyes Open 0.157 0.277 0.894
Vocabulary AF3 Theta Eyes Open 0.042 0.770 0.953
Vocabulary AF3 Alpha Eyes Open 0.176 0.220 0.343
Vocabulary AF3 Low Beta Eyes Open 0.194 0.177 0.394
Vocabulary AF3 Upper Beta Eyes Open 0.218 0.128 0.358
Vocabulary AF3 High Beta Eyes Open 0.130 0.370 0.991
Vocabulary AF3 Gamma Eyes Open 0.080 0.583 0.720
Vocabulary AF4 Delta Eyes Open 0.090 0.536 0.894
Vocabulary AF4 Theta Eyes Open 0.062 0.666 0.953
Vocabulary AF4 Alpha Eyes Open 0.196 0.173 0.343
Vocabulary AF4 Low Beta Eyes Open 0.267 0.061 0.282
Vocabulary AF4 Upper Beta Eyes Open 0.253 0.077 0.269
Vocabulary AF4 High Beta Eyes Open 0.139 0.334 0.991
Vocabulary AF4 Gamma Eyes Open 0.079 0.586 0.720
Vocabulary F3 Delta Eyes Open -0.136 0.346 0.894
Vocabulary F3 Theta Eyes Open 0.042 0.772 0.953
Vocabulary F3 Alpha Eyes Open 0.165 0.253 0.354
Vocabulary F3 Low Beta Eyes Open 0.087 0.549 0.641
Vocabulary F3 Upper Beta Eyes Open 0.182 0.206 0.380
Vocabulary F3 High Beta Eyes Open -0.002 0.991 0.991
Vocabulary F3 Gamma Eyes Open -0.096 0.509 0.720
Vocabulary F4 Delta Eyes Open -0.122 0.398 0.894
Vocabulary F4 Theta Eyes Open 0.043 0.768 0.953
Vocabulary F4 Alpha Eyes Open 0.177 0.218 0.343
Vocabulary F4 Low Beta Eyes Open 0.176 0.222 0.394
Vocabulary F4 Upper Beta Eyes Open 0.201 0.163 0.379
Vocabulary F4 High Beta Eyes Open 0.010 0.946 0.991
Vocabulary F4 Gamma Eyes Open -0.114 0.430 0.720
Vocabulary F7 Delta Eyes Open 0.156 0.279 0.894
Vocabulary F7 Theta Eyes Open 0.087 0.550 0.953
Vocabulary F7 Alpha Eyes Open 0.188 0.192 0.343
Vocabulary F7 Low Beta Eyes Open 0.215 0.133 0.394
Vocabulary F7 Upper Beta Eyes Open 0.258 0.070 0.269
Vocabulary F7 High Beta Eyes Open 0.142 0.325 0.991
Vocabulary F7 Gamma Eyes Open 0.072 0.617 0.720
Vocabulary F8 Delta Eyes Open -0.050 0.729 0.894
Vocabulary F8 Theta Eyes Open 0.001 0.995 0.995
Vocabulary F8 Alpha Eyes Open 0.208 0.148 0.343
Vocabulary F8 Low Beta Eyes Open 0.175 0.225 0.394
Vocabulary F8 Upper Beta Eyes Open 0.178 0.217 0.380
Vocabulary F8 High Beta Eyes Open 0.013 0.928 0.991
Vocabulary F8 Gamma Eyes Open -0.104 0.473 0.720
Vocabulary FC5 Delta Eyes Open -0.073 0.617 0.894
Vocabulary FC5 Theta Eyes Open -0.172 0.233 0.953
Vocabulary FC5 Alpha Eyes Open 0.063 0.664 0.721
Vocabulary FC5 Low Beta Eyes Open 0.049 0.737 0.737
Vocabulary FC5 Upper Beta Eyes Open 0.085 0.559 0.711
Vocabulary FC5 High Beta Eyes Open 0.006 0.970 0.991
Vocabulary FC5 Gamma Eyes Open -0.019 0.897 0.897
Vocabulary FC6 Delta Eyes Open -0.166 0.248 0.894
Vocabulary FC6 Theta Eyes Open -0.037 0.799 0.953
Vocabulary FC6 Alpha Eyes Open 0.110 0.446 0.568
Vocabulary FC6 Low Beta Eyes Open 0.109 0.450 0.610
Vocabulary FC6 Upper Beta Eyes Open 0.054 0.710 0.764
Vocabulary FC6 High Beta Eyes Open -0.018 0.902 0.991
Vocabulary FC6 Gamma Eyes Open -0.101 0.485 0.720
Vocabulary O1 Delta Eyes Open -0.096 0.506 0.894
Vocabulary O1 Theta Eyes Open -0.049 0.734 0.953
Vocabulary O1 Alpha Eyes Open 0.185 0.198 0.343
Vocabulary O1 Low Beta Eyes Open 0.102 0.480 0.610
Vocabulary O1 Upper Beta Eyes Open 0.163 0.258 0.401
Vocabulary O1 High Beta Eyes Open -0.033 0.818 0.991
Vocabulary O1 Gamma Eyes Open -0.108 0.457 0.720
Vocabulary O2 Delta Eyes Open -0.043 0.766 0.894
Vocabulary O2 Theta Eyes Open 0.026 0.859 0.953
Vocabulary O2 Alpha Eyes Open 0.242 0.090 0.343
Vocabulary O2 Low Beta Eyes Open 0.280 0.049 0.282
Vocabulary O2 Upper Beta Eyes Open 0.293 0.039 0.269
Vocabulary O2 High Beta Eyes Open 0.058 0.692 0.991
Vocabulary O2 Gamma Eyes Open -0.081 0.574 0.720
Vocabulary P7 Delta Eyes Open -0.123 0.395 0.894
Vocabulary P7 Theta Eyes Open -0.171 0.234 0.953
Vocabulary P7 Alpha Eyes Open 0.062 0.670 0.721
Vocabulary P7 Low Beta Eyes Open -0.054 0.712 0.737
Vocabulary P7 Upper Beta Eyes Open 0.001 0.994 0.994
Vocabulary P7 High Beta Eyes Open -0.127 0.378 0.991
Vocabulary P7 Gamma Eyes Open -0.155 0.281 0.720
Vocabulary P8 Delta Eyes Open -0.067 0.646 0.894
Vocabulary P8 Theta Eyes Open 0.100 0.491 0.953
Vocabulary P8 Alpha Eyes Open 0.241 0.092 0.343
Vocabulary P8 Low Beta Eyes Open 0.296 0.037 0.282
Vocabulary P8 Upper Beta Eyes Open 0.313 0.027 0.269
Vocabulary P8 High Beta Eyes Open 0.082 0.574 0.991
Vocabulary P8 Gamma Eyes Open -0.049 0.736 0.793
Vocabulary T7 Delta Eyes Open 0.006 0.968 0.968
Vocabulary T7 Theta Eyes Open -0.202 0.160 0.953
Vocabulary T7 Alpha Eyes Open 0.012 0.935 0.935
Vocabulary T7 Low Beta Eyes Open -0.109 0.450 0.610
Vocabulary T7 Upper Beta Eyes Open -0.097 0.502 0.703
Vocabulary T7 High Beta Eyes Open -0.132 0.359 0.991
Vocabulary T7 Gamma Eyes Open -0.108 0.454 0.720
Vocabulary T8 Delta Eyes Open 0.006 0.968 0.968
Vocabulary T8 Theta Eyes Open 0.021 0.885 0.953
Vocabulary T8 Alpha Eyes Open 0.179 0.213 0.343
Vocabulary T8 Low Beta Eyes Open 0.186 0.196 0.394
Vocabulary T8 Upper Beta Eyes Open 0.057 0.693 0.764
Vocabulary T8 High Beta Eyes Open -0.098 0.497 0.991
Vocabulary T8 Gamma Eyes Open -0.103 0.477 0.720

And here is the corresponding distributions of r and p values. The dashed lines correspond to a significant threshold of p < .05 on either a r or a p-value scale.

ggplot(Rdata_eo, aes(x = Band, y = r, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", 
               alpha=0.5, geom = "errorbar") +
  facet_wrap(~ material) +
  ggtitle("Correlation, by Band") +
  ylab("r value") +
  annotate("segment", x=-Inf, xend=Inf, y=0.28, yend = 0.28, lty=2) +
  annotate("segment", x=-Inf, xend=Inf, y=-0.28, yend = -0.28, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  theme_pander()

ggplot(Rdata_eo, aes(x = Band, y = p, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", 
               alpha=0.5, geom = "errorbar") +
  ggtitle("p-value, by Band") +
  ylab("p value") +
  facet_wrap(~ material) +
  annotate("segment", x=-Inf, xend=Inf, y=0.05, yend = 0.05, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  scale_y_log10() +
  theme_pander()

Now, we might wonder how much the correlations change from eyes closed to eyes open recordings. The best way to chack is to visualize then as a scatterplot:

Rdata <- bind_rows(Rdata_ec, Rdata_eo)
wrdata <- pivot_wider(Rdata, id_cols = c("material", "Channel", "Band"),
                      names_from = Recording, 
                      values_from = r)

ggplot(wrdata, aes(x = `Eyes Open`, y = `Eyes Closed`, 
                   col = material)) +
  annotate("segment", x = 0, y=-Inf, xend=0, 
           yend=Inf, col="grey", lwd=1, lty=1) +
  annotate("segment", x = -Inf, y=0, xend=Inf, 
           yend=0, col="grey", lwd=1, lty=1) +
  
  geom_point(alpha=0.5, size=3) +
  annotate("segment", x = -0.25, y=-0.25, xend=0.5, 
           yend=0.5, col="black", lwd=1, lty=2) +
  theme(plot.title = element_text(hjust = 0.5)) +
  #geom_smooth(method = "lm", formula = y ~ x, 
  #            fullrange = T, lwd = 1) +
  scale_color_brewer(palette="Dark2") +
  ggtitle("Correlation Coefficients Across Recordings") +
  theme_pander()

Let’s tabulate only those channels that survive the FDR correction:

survivors <- Rdata_eo %>%
  filter(q < 0.05)
  
survivors %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
material Channel Band Recording r p q
Maps AF3 Upper Beta Eyes Open 0.383 0.006 0.021
Maps AF4 Upper Beta Eyes Open 0.383 0.006 0.021
Maps F3 Upper Beta Eyes Open 0.326 0.021 0.037
Maps F4 Upper Beta Eyes Open 0.338 0.016 0.034
Maps F7 Upper Beta Eyes Open 0.428 0.002 0.013
Maps F8 Upper Beta Eyes Open 0.336 0.017 0.034
Maps O2 Upper Beta Eyes Open 0.344 0.015 0.034
Maps P8 Upper Beta Eyes Open 0.436 0.002 0.013

In the case of eyes-open data, we have a much larger set of channels that survive correction for multiple comparison. Although these correlations generally reflect the same trends we have seen the eyes-closed data (i.e., significant positive correlations for prefrontal and right parietal sites), they are now stroger for visual (not verbal) materials and centered at at slightly higher frequency (upper beta instead of low beta).

Here is an overview of the scatterplots and correlations:

focus <- Fdata_eo %>%
  filter(material == "Maps", 
         Channel %in% survivors$Channel,
         Band == "Upper Beta") %>%
  rename(Alpha = alpha)

ggplot(focus, aes(x=Alpha, y=Power)) +
  geom_point(size=4, alpha=0.5, col="black") + 
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() +   
  scale_y_continuous() + 
  ggtitle("Rate of Forgetting and Eyes-Open Beta Power") +
  xlab(expression(paste(alpha, " Maps"))) +
  ylab("Upper Beta Power (15-18 Hz)") +
  facet_wrap(~ Channel, scales="free_y") +
  geom_text(data=survivors, col="red",
            mapping=aes(x=0.25, y=9.5, 
                        label= paste("r =", round(r, 2)))) +
  theme(panel.spacing = unit(1.5, "lines")) +
  theme(plot.title = element_text(hjust = 0.5))

Predicting Rate of Forgetting with Group Lasso

So far, we have shown that neurophysiological signals exists that are correlated with the rate of forgetting. This, however, does not directly answer wheter we could “read” an individual’s rate of forgetting from their EEG data directy.

Group LASSO

To do so, we need to use a ML approach. Specifically, here we will use group LASSO regresso. Group LASSO is a way to combine regression, validation, and feature selection into a single approach.

LASSO works like a normal regression, with the additional constraint that the \(\beta\) regressors need to minimize the sum of the residual sum of squares (RSS) and the quantity \(- \lambda||\beta||_{1}\), where \(||\beta||_{1}\) is the first-order norm (the sum of absolute values of each regressor weight). For \(\lambda = 0\), LASSO reduces to normal regression but, as the value of \(\lambda\) grows, more and more regressors are “pushed” to zero and removed from the pool. In group LASSO, regressors are combined and selected by groups, which is a natural choice when natural groupings exists, as in the cases of frequency bands.

As a demo, we will pilot this approach to eyes-closed recordings and the combined rates of forgetting.an

Data Preparation

First, we need to transform the data appropriately, converting from the long to the wide format, which is better suited for regression:

library(gglasso)
wfdata <- Fdata_ec %>% 
  pivot_wider(id_cols = c("Subject", "Channel", 
                          "Band", "material", "alpha"), 
              names_from = c("Channel", "Band"), 
              values_from = c("Power"))

wide_data_ec <- wfdata %>%
  pivot_wider(names_from = c("material"), 
              values_from = "alpha")

From the wide format data, we are going to define the data that will enter the regression, the matrix of observations X and the vector of variables Y. The matrix of observations X is 50 x 98 matrix, corresponding to 50 participants x 98 possible power values (14 channels x 7 frequency bands). The Y vector is made of 50 observed rates of forgetting. To make the prediction reliable, we will average the rates of forgetting across materials (Maps + Vocabulary)

channel_names <- c("AF3", "AF4", "F3", "F4", 
                   "F7", "F8", "FC5", "FC6", 
                   "T7", "T8", "P7", "P8", 
                   "O1", "O2")


band_names <- c("Delta", "Theta", "Alpha", "Low Beta", 
                "Upper Beta", "High Beta", "Gamma")

by_band_order <- as.vector(
  outer(channel_names, 
        band_names, 
        function(x, y){
          paste(x, y, sep="_")
          }))

# Predictors
X <- wide_data_ec[,c(-1, -100, -101)] %>%
  #select(ends_with("Beta")) %>%
  as.matrix()

X <- X[,by_band_order]
# We need to reorder X by frequency bands


# Maps
Ymaps <- wide_data_ec[,c(101)] %>%
  as.matrix()

# Vocab
Yvoc <- wide_data_ec[,c(101)] %>%
  as.matrix()

# Average across materials to obrtain stable estimates
Y <- (Ymaps + Yvoc)/2

# Let's save the mean (will be halpful later)
uY <- mean(Y)

# Center Y to remove dominant effect of intercept
Y <- Y - uY

Finally, we need to define our set of grouping variables, defined as a vector of numbers that identify each group. We are going to create two such vectors, one for the freqency bands and one for the channels. In this example we are going to use bands as the grouping variables.

bands <- rep(1:7, each=14)
channels <- sort(rep(1:14, 7))

Fitting the LASSO

Now, with that in place, we can fit the LASSO model. The LASSO model is fit much like a normal regression, with an additional partameter \(\lambda\).

Because we do not know a priori which value of \(lambda\) to pick, the canonical solution is to fit fit different models for different values of \(\lambda\).

fit <- gglasso(x = X, y = Y, 
               group = bands,
               loss = 'ls', 
               nlambda=100,
               lambda.factor = 0.001)
coef.mat=fit$beta

Finding the optimal \(\lambda\) with Cross-Validation

Our goal, of course, is generalizability across participants. Although the fit model is important, what we really care is a different measure, that of cross-validation. For cross-validation, we are going to use the n-fold approach, in which the data is randomly split so that 1/n of it is saved for test. (if n is identical to the number of participants, it becomes the “Leave-One-Out” validation scheme). Here, we used the Leave-One-Out method, which guarantees that the results are the same every time this script is ruin (because there only one possible permutation of the data). So, n = 50.

fit.cv <- cv.gglasso(x = X,
                     y = Y,
                     group = bands,
                     pred.loss = "L1",
                     nlambda = 100,
                     lambda.factor = 0.001,
                     nfolds = 50)

Visualizing Validation Results

Here are the fit and cross-validation plots, with the vertical lines marking the points corresponding to \(\lambda_{min}\) (optimal for cross-validation) and
\(\lambda_{1SE}\) (more conservative model). Because we have used the Leave-One-Out validation scheme, we are can use the the more aggressive \(\lambda_{min}\) estimate, since our sample of the folds was exhaustive.

# Best lambda

lmin <- fit.cv$lambda.min  # Optimal lambda, maximizes generalization
l1se <- fit.cv$lambda.1se  # Conservative model, min + 1SE

plot(fit.cv, 
     pch=20,
     lwd=2)

plot(fit, group=T, 
     col = inferno(length(unique(fit$group)) + 1), 
     lwd = 2)

legend(x="topright", 
       legend=c("Delta", "Theta", "Alpha",
                "Low Beta", "Upper Beta", "High Beta",
                "Gamma"),
       bty="n",
       ncol=2,
       col=inferno(length(unique(fit$group)) + 1),
       lty=1, lwd=2)
abline(v=log(lmin), lty=2)
abline(v=log(l1se), lty=2)

Using the Lambda min value, we can now compare predicted vs. observed values.

prediction <- predict(object = fit,
                    newx = X,
                    s = lmin,
                    type='link')

observed <- data.frame(Subject = wide_data_ec$Subject,
                       RoF = as_vector(Y) + uY,
                       Condition = "Observed") %>%
  rename(RoF = Maps)

predicted <- data.frame(Subject = wide_data_ec$Subject,
                       RoF = as_vector(prediction) + uY,
                       Condition = "Predicted") %>%
  rename(RoF = X1)

comparison <- as_tibble(rbind(observed, predicted))

plt = cbind(Y,
            predict(object = fit,
                    newx = X,
                    s = lmin,
                    type='link'))

# matplot(plt,
#         main="Predicted vs Actual",
#         type='b',
#         lwd=2,
#         pch="o",
#         col=viridis(3),
#         xlab="Subject",
#         ylab="Rate of Forgetting (alpha)")

ggplot(comparison, aes(x=Subject, y=RoF), col=Condition) +
  geom_point(aes(col=Condition)) +
  geom_line(alpha=0.5, lty=2) +
  annotate("text", x= 40, y = 0.20, 
           label=(paste("r =", 
                        round(cor(Y, prediction),
                              2)))) +
  theme_pander() +
  ylab("Rate of Forgetting")+
  ggtitle("Cross-Validation:\nPredicted vs. Observed Rates of Forgetting") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) 

wcomparison <- comparison %>%
  pivot_wider(names_from=Condition, values_from = RoF)

p <- ggplot(wcomparison, aes(x = Predicted, y = Observed)) +
  geom_point(alpha=0.5, size=4) +
  geom_abline(slope = 1, intercept = 0,
              lwd=1, col="blue", lty=2) +
  annotate("text", x= 0.40, y = 0.25, 
           label=(paste("r =", 
                        round(cor(Y, prediction),
                              2)))) +
  theme_pander() +
  #coord_fixed() +
  coord_equal(ratio=1) +
  ylab("Rate of Forgetting (Observed)") +
  xlab("Rate of Forgetting (Predicted)") +
  ggtitle("Cross-Validation") 

ggMarginal(p, col="white", fill="black", alpha=0.5,
           type="hist", bins=13)

In terms of overall error, the root mean squared error between predicted and observed Rate of Forgetting is RMSE = 0.024

And now, let’s print the full table of beta regression weights (for inspectio and analysis) at the desiored value of \(\lambda_{min}\). We can also save the weights for future analysis.

get_channel <- function(s) {
  str_split(s, "_", 2)[[1]][1]
}

get_band <- function(s) {
  str_split(s, "_", 2)[[1]][2]
}

vget_channel <- Vectorize(get_channel)
vget_band <- Vectorize(get_band)


coef(fit.cv, 
     s=lmin) -> betas
  
bdf <- data.frame("Predictor" = rownames(betas),
                  "Beta" = as.numeric(betas)) %>%
  mutate(Channel = vget_channel(Predictor),
         Band = vget_band(Predictor)) 

bdf %>%
  filter(Beta != 0) %>%
  xtable() %>%
  kable(digits = 5) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Predictor Beta Channel Band
(Intercept) 0.00955 (Intercept) NA
AF3_Delta -0.00659 AF3 Delta
AF4_Delta 0.00294 AF4 Delta
F3_Delta -0.02081 F3 Delta
F4_Delta 0.03104 F4 Delta
F7_Delta 0.01171 F7 Delta
F8_Delta -0.01349 F8 Delta
FC5_Delta -0.01473 FC5 Delta
FC6_Delta -0.01259 FC6 Delta
T7_Delta -0.00654 T7 Delta
T8_Delta 0.01771 T8 Delta
P7_Delta 0.00053 P7 Delta
P8_Delta 0.02433 P8 Delta
O1_Delta 0.00076 O1 Delta
O2_Delta 0.01153 O2 Delta
AF3_Theta -0.00014 AF3 Theta
AF4_Theta -0.00003 AF4 Theta
F3_Theta -0.00103 F3 Theta
F4_Theta -0.00013 F4 Theta
F7_Theta -0.00101 F7 Theta
F8_Theta -0.00092 F8 Theta
FC5_Theta -0.00081 FC5 Theta
FC6_Theta 0.00035 FC6 Theta
T7_Theta -0.00324 T7 Theta
T8_Theta 0.00128 T8 Theta
P7_Theta -0.00182 P7 Theta
P8_Theta 0.00146 P8 Theta
O1_Theta 0.00001 O1 Theta
O2_Theta 0.00095 O2 Theta
AF3_Alpha -0.00175 AF3 Alpha
AF4_Alpha 0.00352 AF4 Alpha
F3_Alpha 0.00034 F3 Alpha
F4_Alpha 0.00289 F4 Alpha
F7_Alpha 0.00115 F7 Alpha
F8_Alpha 0.00520 F8 Alpha
FC5_Alpha 0.00101 FC5 Alpha
FC6_Alpha 0.01122 FC6 Alpha
T7_Alpha -0.00497 T7 Alpha
T8_Alpha 0.01223 T8 Alpha
P7_Alpha -0.01927 P7 Alpha
P8_Alpha -0.00379 P8 Alpha
O1_Alpha -0.00753 O1 Alpha
O2_Alpha 0.00392 O2 Alpha
AF3_Upper Beta 0.00068 AF3 Upper Beta
AF4_Upper Beta 0.00089 AF4 Upper Beta
F3_Upper Beta 0.00009 F3 Upper Beta
F4_Upper Beta 0.00030 F4 Upper Beta
F7_Upper Beta 0.00042 F7 Upper Beta
F8_Upper Beta 0.00064 F8 Upper Beta
FC5_Upper Beta 0.00014 FC5 Upper Beta
FC6_Upper Beta 0.00069 FC6 Upper Beta
T7_Upper Beta 0.00046 T7 Upper Beta
T8_Upper Beta 0.00027 T8 Upper Beta
P7_Upper Beta 0.00023 P7 Upper Beta
P8_Upper Beta 0.00002 P8 Upper Beta
O1_Upper Beta 0.00058 O1 Upper Beta
O2_Upper Beta 0.00066 O2 Upper Beta
AF3_Gamma 0.00144 AF3 Gamma
AF4_Gamma 0.00539 AF4 Gamma
F3_Gamma -0.01068 F3 Gamma
F4_Gamma -0.00540 F4 Gamma
F7_Gamma -0.01039 F7 Gamma
F8_Gamma -0.00435 F8 Gamma
FC5_Gamma -0.00148 FC5 Gamma
FC6_Gamma -0.00273 FC6 Gamma
T7_Gamma 0.00736 T7 Gamma
T8_Gamma -0.01448 T8 Gamma
P7_Gamma -0.00609 P7 Gamma
P8_Gamma -0.00418 P8 Gamma
O1_Gamma 0.00834 O1 Gamma
O2_Gamma 0.01031 O2 Gamma
write_csv(bdf, "group_lasso_weights.csv", col_names = T)

And here is the visualization of the corresponding weights as a topomap:

Topomap of Beta weights